% Draw coefficients from Multivariate Regressions models
%  see e.g. Zellner (1971), pp. 224-240
% See equation 8.14 on page 227 for posterior with diffuse priors and
% equations 8.71 - 8.73 for posterior with informative priors

function [b,varargout]=multiReg(Y,X,Sigma,V_prior,B_prior,sv)
n=size(Y,2);

if sv>1
    beta_hat = vec((X'*X)^(-1)*X'*Y);
    
    B = eye(size(V_prior))/( eye(size(V_prior))/V_prior + kron((Sigma)^(-1),X'*X))*( eye(size(V_prior))/V_prior*B_prior + kron((Sigma)^(-1),X'*X)*beta_hat);
    V = eye(size(V_prior))/(eye(size(V_prior))/V_prior+kron((Sigma)^(-1),X'*X));
    
else
    X = kron(eye(n),X);
    Omegat = gen_omegat(Sigma,n);
    
    beta_hat = vec((X'/Omegat*X)^(-1)*X'/Omegat*Y);
    
    B = eye(size(V_prior))/( eye(size(V_prior))/V_prior + X'/Omegat*X )*( eye(size(V_prior))/V_prior*B_prior + X'/Omegat*X*beta_hat);
    V = eye(size(V_prior))/(eye(size(V_prior))/V_prior+ X'/Omegat*X);
    
end
b = mvnrnd1(B,V,1);

varargout(1) = {B};
varargout(2) = {V};


